An exact formalism for quench dynamics 
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We describe a formulation for studying the quench dynamics of integrable systems generalizing an 
approach by Yudson. We study the evolution of the Lieb-Liniger model, a gas of interacting bosons 
moving on the continuous infinite line and interacting via a short range potential. The formalism 
allows us to quench the system from any initial state. We find that for any value of repulsive 
coupling independently of the initial state the system asymptotes towards a strongly repulsive gas, 
while for any value of attractive coupling, the system forms a maximal bound state that dominates 
at longer times. In either case the system equilibrates but does not thermalize. We compare this to 
quenches in a Bose-Hubbard lattice and show that there, initial states determine long-time dynamics 
independent of the sign of the coupling. 



I. INTRODUCTION 

Noncquilibrium processes occur in fields as diverse as 
metallurgy and cell biology - in fact, most physical pro- 
cesses that occur in nature are dynamical. It is there- 
fore important to understand nonequilibrium phenomena 
from a variety of perspectives and in several different con- 
texts. 

In the context of physics, the dynamics of nonequi- 
hbriuni processes has been studied since the early days 
of thermodynamics (see Ref. [l| for a historical introduc- 
tion to the subject). More recently, the study of quan- 
tum nonequilibrium processes has received a boost from 
highly tunable experimental systems, in particular ultra- 
cold atomic gases, that can be well isolated from their 
environment. They provide a testing ground for theo- 
ries of quantum nonequilibrium behavior and have led to 
observations of interesting new phenomena ^ . 

Unlike thermodynamics, there exists no general frame- 
work to understand far from equilibrium behavior. Dif- 
ferent problems need different approaches and in most 
cases we can only understand the behavior of physical 
observables in limited regimes of parameters or windows 
of time. However, a good sense of the complexity of the 
problem is emerging from a multitude of studies on sev- 
eral types of systems analytically, computationally, and 
exp er iment ally. 

Within the domain of quantum phenomena, transport 
and quenches are the primary means of studying nonequi- 
librium physics. Transport phenomena include for exam- 
ple, transient and steady-state currents in devices, study 
of quantum impurity models and their conduction prop- 
erties, quantum hall edge states, and the surface states 
of topological insulators. Quenches, on the other hand, 
allow us to study phenomena like thermalization and re- 
laxation of physical observables by strongly disturbing an 
equilibrium system and watching it evolve Q • Cold atom 
systems, as noted, are particularly amenable to quenches 
- we briefly describe these in section |lll On the theo- 
retical side, with advances in computational techniques, 
it is possible to calculate the time evolution of various 



observables of larger systems and study the underlying 
physics 0-01 • However, these techniques are not suitable 
for very large systems or continuum models. 

In this article, we are interested in studying the quench 
dynamics of systems described by integrable models. 
These models possess an infinite number of conserved 
quantities, and this property is expected to reveal itself 
in their dynamics [J, la, |8| . A generic model would access 
all of the constant energy surface in phase space as it 
evolves, but this may not be the case for an integrable 
model. 

The eigenstates of integrable models can be obtained 
via the Bethe Ansatz technique [9[. The Bethe Ansatz 
has proved to be of tremendous use in studying the 
ground states and thermodynamics of such models |10l - 
Il3| . However, in spite of knowing the eigenstates of the 
Hamiltonian the dynamics still remains a complicated 
problem. In the remainder of this paper, we elucidate a 
framework for dynamics introduced by Yudson [14l Il5l | . 
We introduce some generalizations and use it to under- 
stand the quench dynamics of a gas of bosons with at- 
tractive or repulsive short ran ge i nteractions in the con- 
text of the Lieb-Liniger model [Ig . This approach allows 
analytical calculations that are essential for a complete 
understanding of the out-of-equilibrium properties of the 
system. 

The article is a follow up to Ref. (iT). Some of the 
formulas are derived in detail, and some of the plots are 
repeated for completeness. We also correct an error in 
the earlier work. 

Before getting into details we discuss some general is- 
sues concerning the quench dynamics of isolated systems 
defined on the infinite line. 



A. Dynamics of Id isolated many-body systems 

While studying the thermodynamic properties of a 
quantum system, one needs to enumerate and classify 
the eigenstates of the Hamiltonian in order to construct 
the partition function. To achieve this, some finite vol- 



umc boundary conditions (BC) arc typically imposed — 
either periodic BC to maintain translation invariance or 
open BC when the system has physical ends. One can 
then identify the ground state and the low lying excita- 
tions that dominate the low-temperature physics. In the 
thermodynamic limit "^^^^ = p, the limit of very large 
systems with large number of particles at finite density, 
the effect of the boundary condition is negligible and we 
expect the results to be universally valid. 

When the system is out of equilibrium, a different set 
of issues arises. We shall consider here the process of a 
"quantum quench" where one studies the time evolution 
of a system after a sudden change in the parameters of 
the Hamiltonian governing the system. To be precise, 
one assumes that the system starts in some stationary 
state l^o)- This stationary state can be thought of as 
the ground state of some (interacting) Hamiltonian Hq. 
Following the quench at i = 0, the system evolves in time 
under the influence of a new Hamiltonian H which may 
differ from Hq in many ways. One may add an interac- 
tion, change an interaction coupling constant, apply or 
remove an external potential or increase the size of the 
system. Further, the quench can be sudden, i.e., over a 
time window much shorter than other time scales in the 
system, driven at a constant rate or with a time depen- 
dent ramp. 

In this paper we shall concentrate on sudden quenches 
where the initial state I^Pq) describes a system in a finite 
region of space with a particular density profile: lattice- 
like, or condensate-like (see Fig. [T]). Under the effect of 
the quenched Hamiltonian the system evolves as j^Oi t) ~ 




FIG. 1. Initial states, (a) For ^ ^ 1, we have a lattice like 
state, l^iatt). (b) For a = 0, we have a condensate like state 
l^cond), o" determines the spread. 

To compute the evolution, it is convenient to expand 
the initial state in the eigenbasis of the evolution Hamil- 
tonian, 



where \n) are the eigenstates of H and C„ = {n\'^o) 
are the overlaps with the initial state, determining the 
weights with which different eigenstates contribute to the 
time evolution: 

|vI/o,t)-^e— "*C„|n). (L2) 

{"} 

The evolution of observables is then given by, 

(d(t))*„ = (vl/o,i|0|vl/o,i) 

= ^ e-^('"-^")*C:C„(ri|d|m), (1-3) 

{n,m} 

with observables that may be local operators, correlation 
functions, currents or global quantities such as entangle- 
ments. 

The time evolution is characterized by the energy of 
the initial state. 
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which is conserved throughout the evolution, specifying 
the energy surface on which the system moves. This sur- 
face is determined by the initial state through the over- 
laps C„. Unlike the situation in thermodynamics where 
the ground state and low-lying excitations play a central 
role, this is not the case out-of-equilibrium. A quench 
puts energy into the system which the isolated system 
cannot dissipate and it cannot relax to its ground state. 
Rather, the eigenstates that contribute to the dynamics 
depend strongly on the initial state via the overlaps C„ 
(see Fig. HI). 

A vivid illustration comes from comparing quenches 
in systems that differ in the sign of the interaction. In 
the Bose-Hubbard model and the XXZ model it has been 
observed that the sign of the interaction plays no role in 
the quench dynamics [l^, [l^, even though the ground 
states that correspond to the different signs are very dif- 
ferent. For example, for the XXZ magnet, the ground 
state is either ferromagnetic or Ncel ordered (RVB in Id) 
depending on the sign of the anisotropy A. In Appendix 
[M we show results for the Bose-Hubbard model, and pro- 
vide an argument for this observation. The Lieb-Liniger 
model, whose quench dynamics we describe here, on the 
other hand shows very different behavior, reaching an 
long time equilibrium state that depends mainly on the 
sign of the interaction. 

In the experiments that we seek to describe, a system 
of N bosons is initially confined to a region of space of 
size L and then allowed to evolve on the infinite line while 
interacting with short range interactions. It is important 
to first understand the time-scales of the phenomena that 
we are studying. There are two main types of time scales 
here. One is determined by the initial condition (spa- 
tial extent, overlap of nearby wave- functions) , and the 
other by the parameters of the quenched system (mass, 
interaction strength). 
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FIG. 2. Difference between quencli dynamics and thermody- 
namics. After a quench, the system probes high energy states 
and does not necessarily relax to the ground state. In ther- 
modynamics, we minimize the energy (or free energy) of a 
system and probe the region near the ground state. 



For an extended system where we start with a locally 
uniform density (sec Fig.lTJi), wc expect the dynamics to 
be in the constant density regime as long as t ^ —, v be- 
ing the characteristic velocity of propagation. Although 
the low energy thermodynamics of a constant density 
Bose gas can be described by a Luttinger liquid [20|, we 
expect the collective excitations of the quenched system 
to behave as a highly excited Liquid since the initial state 
is far from the ground state. It is also possible that de- 
pending on the energy density equcnch/^j the Luttinger 
liquid description may break down altogether. 

The other time scale that enters the description of 
nonequilibrium dynamics is the interaction time scale, r, 
a measure of the time it takes the interactions to develop 



fully: 



for the Lieb-Liniger model [21| . Assum- 



ing L is large enough so that r ^ — , we expect a fully 
interacting regime to be operative at times beyond the 
interactions scale until t ^ — and the density of the sys- 
tem can no longer be considered constant, diminishing 
with time as the system expands. In the Lieb-Liniger 
model, this leads to an effective increase in the coupling 
constant which manifests itself as fermionization for re- 
pulsive interaction and bound-state correlations in the 
case of attractive interactions. Thus the main operation 
of the interaction occurs in the time range r < i < — , 
over which the wave function rearranges and after which 
the system is dilute and freely expands. For the case 
L = oo, free expansion is not present. Figure [3] sum- 
marizes the different time-scales involved in a dynamical 
situation. 

We shall also consider initial conditions where the 
bosons are "condensed in space", occupying the same 
single particle state characterized by some scale a (see 
Fig. [H)). In this case the short time dynamics is not 
present, the time scales at which we can measure the 
system are typically much larger than -^ and we expect 



excited Luttinger liquid 



free expansion 



' t ^ 

fully interacting regime 



FIG. 3. Time scales involved in quench dynamics, r is an 
intrinsic time scale that depends on the interaction strength. 
L/v is a characteristic time at which the system sees the finite 
extent. 



the dynamics to be in the strongly interacting and ex- 
panding regime. 



B. Quench dynamics and the Bethe Ansatz 

To carry out the computation of the quench dynam- 
ics we need to know the eigenstates of the propagating 
Hamiltonian. The Bethe Ansatz approach is helpful in 
this respect as it provides us with the eigenstates of a 
large class of interacting one dimensional Hamiltonians. 
Many of the Hamiltonians that can be thus be solved are 
of fundamental importance in condensed matter physics 
and have been proposed to study various experimental 
situations. A partial list includes the Heisenberg chain 
(magnetism), the Hubbard model (strong-correlations), 
the Lieb-Liniger model (cold atoms in optical traps), 
the Kondo model, and the Anderson model (impuri- 
ties in metals, quantum dots) 0, [lll - [l3l \WL |22| . For a 
Hamiltonian to possess eigenstates that are given in the 
form of a Bethe Ansatz it must have the property that 
multi-particle interactions can be consistently factorized 
into series of two particle interactions, all of them being 
equivalent Q. 

While originally formulated to understand the Heisen- 
berg chain Q, the technology has been studied exten- 
sively and recast into a more sophisticated algebraic for- 
mulation 12J, [23 ■ The usual focus of the Bethe Ansatz 
approach has been on the thermodynamic properties of 
the system: determining the spectrum, the free energy, 
and susceptibilities. Also, considerable efforts were made 
to compute correlation functions |25l. |26|. In particular, 
the Algebraic Bethe Ansatz has been used in conjunc- 
tion with numerical methods to calculate so-called form 
factors which allows access to the dynamical structure 
functions d 0,113 ■ 

Using the Bethe Ansatz to extract dynamics one en- 
counters a triply complicated problem — the first being 
to obtain the full spectrum of the Hamiltonian, the sec- 
ond to calculate overlaps, and the third to carry out the 
sum, which in some cases involves sums over large sets 
of different "configurations" of states. The overlaps are 
particularly difficult to evaluate due to the complicated 
nature of the the Bethe eigenstates and their normaliza- 



tion. The problem is more pronounced in the far from 
equihbrium ease of a quench when the state we start with 
suddenly finds itself far away from the eigenstates of the 
new Hamiltonian [see eq. p.2p ]. and all the eigenstates 
have non-trivial weights in the time-evolution. In all but 
the simplest cases, the problem is non-perturbative and 
the existing analytical techniques are not suited for a di- 
rect application to such a situation. 

Often quench calculations are carried out in finite vol- 
ume and the infinite volume limit is taken at the end of 
the calculation. This may be necessary for thermody- 
namic calculations, as noted, but not for quench dynam- 
ics (see sec. II A[) . We shall instead carry out the quench 
directly in the infinite volume limit in which case the 
Yudson representation allows us to carry out the calcu- 
lations in an efficient way, doing away with some of the 
difficulties mentioned above by not requiring any infor- 
mation about the spectrum, and by using integration as 
opposed to discrete summation. 

More explicitly, the Schrodinger equation for N bosons 
H\X) = e(A)|A) is satisfied for any value of the mo- 
menta {Xj,j = !,■•• , N} if no boundary conditions 
are imposed. The initial state can then be written as 
l^o) = J d^ ^ C\ I A), with the integration over A re- 
placing the summation. This is akin to summing over 
an over-complete basis, the relevant elements in the sum 
being automatically picked up by the overlap with the 
physical initial state. 

It is important to note, however, that while the spec- 
trum of an infinite system is continuous, it can still be 
very complicated. This is indeed the case with several in- 
tegrable models, including the Lieb-Liniger model, where 
the analytic structure of the 5-matrix (the momentum 
dependent phase picked up when two particles cross) 
may permit momenta in the form of complex-conjugate 
pairs signifying bound states. In the formalism we em- 
ploy, such states are taken into account by appropriately 
choosing the contours of integration in the complex plane. 
This procedure is described in detail. 

The remainder of this article is organized as follows. 
In section |TT1 we briefly discuss experiments with cold 
atoms. In sections IlIII and ITVl we will describe the Lieb- 
Liniger model and Yudson representation for this model, 
and show how we can use it to calculate the time evolu- 
tion of an arbitrary initial state. Section |V] studies the 
two particle case in detail. We then go on to calculate the 
evolution of the noise correlations for both the repulsive 
and the attractive gas at long times in sec. I VI I We con- 
clude with a conjecture about the nature of equilibration 
and thermalization and end with a description of ongo- 
ing work and future directions in sec. I VIII Appendix |X] 
discusses a quench in the Bose-Hubbard model, and why 
the sign of the interaction doesn't affect the quench dy- 
namics. 



II. EXPERIMENTS WITH COLD ATOMS 

The experimental arena to which our calculations ap- 
ply are ultracold atomic gases in laser traps, which have 
become a powerful system for exploring nonequilibrium 
phenomena [^, l28| . The systems are formed by trapping 
a gas of atoms using standing light waves made by lasers. 
The gases are cooled evaporatively and are well isolated 
from any thermal baths making them ideal for study- 
ing relaxation and thermalization in isolated quantum 
systems. The interactions between the particles, the po- 
tentials, and their statistics can be controlled by the use 
of external magnetic and electric fields, tuning the op- 
tical lattice, and loading different atoms into the traps. 
Systems with mobile impurities can also be studied by 
loading two or more different species of atoms into the 
lattices. Lattices can be three dimensional or can be 
made quasi-ld or 2d by using confining potentials. The 
typical relaxation and evaporation time scales in these 
systems are in the milliseconds. This makes measure- 
ment easier than in solid state systems. It also allows for 
sudden quenches. Disorder is also largely absent, unless 
introduced. 

Tuning the parameters allows the study of superfluid 
behavior, Mott insulators, spin chains and so on. Such 
a gas trapped by lasers and cooled to nano-Kelvin tem- 
peratures can be quenched by suddenly changing the in- 
teraction between the molecules and the external trap- 
ping potential. Evolution can be globally observed by 
imaging the gas, and the time evolution of densities and 
correlation functions can be obtained from these im- 
ages gSlli,!!^. 

In one dimension, which is of particular interested to 
us, the typical models that are used to study these sys- 
tems are the Bose-Hubbard model, the XXZ model, the 
Sine-Gordon model and the Lieb-Liniger model. Each 
of these models studies a different regime of the gas. 
The Bose-Hubbard model is optimal for atoms hopping 
on a one dimensional lattice. A particular limit of the 
Bose-Hubbard model can be mapped to the XXZ spin 
chain |19| which is intcgrable. The continuum gas is cap- 
tured by the Lieb-Liniger model. 

In this article, we will study the long time dynamics 
of the Lieb-Liniger model, also an integrable model as 
mentioned. In this context, it is an important question 
as to how "integrable" a particular experimental real- 
ization is. Often the experimental setup maintains an 
external trapping potential and including such a poten- 
tial in an integrable Hamiltonian may render the system 
non-integrable. Experimentally, such potentials need to 
be eliminated to the extent possible. This can be par- 
tially achieved by using blue/red detuned lasers to create 
a flat potential well. As has been shown in Ref. [^, the 
dynamics in a particular experiment very closely resem- 
bles what we expect from an integrable model, and it is 
believed that we can indeed create integrable systems to 
a close approximation. This also opens up the question 
of how far from integrability do we need to be in order to 



see the effects of integrability breaking. We discuss this 
point in some detail in the conclusions. 



III. THE LIEB-LINIGER MODEL 

Bosons in one dimensional traps interact via short 
range potentials, which can be well approximated by a 
(5-function interaction. This model was solved in 1963 
by Lieb and Liniger [16[ who originally introduced it to 
overcome shortcomings of other models and approaches 
for understanding quantum gases and liquids, and to pro- 
vide a rigorous result to test perturbation theory against. 
In particular, they sought to improve upon Girardeau's 
hard-core boson model [3l| by providing a tunable pa- 
rameter and better model a low density gas, perhaps ex- 
tensible to higher dimensions. The Schrodinger equation 
for the model is also commonly known as the Non-linear 
Schrodinger equation and has been extensively studied 
both classically and quantum mechanically. 

The Hamiltonian is given by 



H 



[db'' {x)db{x) +cb\x)b{xp{x)b{x)], (III.l) 



where 6 is a bosonic field and c is the interaction strength. 
The mass has been set to 1/2. The action for the model 
~ J^ f[dt — d'^]. Time therefore has the dimension of 
(length)^. The coupling constant c has the dimensions of 
length. 

The model is integrable and the eigenstates take the 
Bethe Ansatz form, 

|A) = NiX) [HZf^iX, - \,)l[e^'^^'^b\x,m, 

(III.2) 



i<3 



where N{X) is a normalization factor determined by a 
particular solution, and 



ZUz) 



Z — iCa^iXi — Xj) 



(III.3) 



incorporates the two particle ^-matrix, Sij = ^'_^^_^^ 
occurring when two bosons cross. The above state satis- 
fies 



iJ|A)=^A2|A) 



(III.4) 



for any value of the momenta A, which, depending on 
the sign of c, may be pure real or form complex pairs. In 
our work, we study the evolution dynamics of the model 
on the infinite line and need not solve for explicit distri- 
butions of the A that characterize the low lying energy 
eigenstates, as discussed in section \LK\ 



IV. YUDSON REPRESENTATION 

In 1985, V. I. Yudson presented a new approach to 
time evolve the Dicke model (a model for superradiance 
in quantum optics [32|) considered on an infinite line |14| . 
The dynamics in certain cases was extracted in closed 
form with much less work than previously required, and 
in some cases where it was even impossible with earlier 
methods. The core of the method is to bypass the labori- 
ous sum over momenta using an appropriately chosen set 
of contours and integrating over momentum variables in 
the complex plane. It is applicable in its original form to 
models with a particular pole structure in the two parti- 
cle S'-matrix, and a linear spectrum. We generalize the 
approach to the case of the quadratic spectrum and apply 
it to the study of quantum quenches. 

As discussed earlier, in order to carry out the quench 
of a system given at i = in a state jVPo) one naturally 
proceeds by introducing a "unity" in terms of a complete 
set of eigenstates and then apply the evolution operator- 



1^0, t) 



-iHt Y^ 



^|A>(A|*o) 
{>■} 



{A} 



e-«(^)*|A)(A|vl/o) 



(IV.l) 



The Yudson representation overcomes the difficulties in 
carrying out this sum by using an integral representa- 
tion for the complete basis directly in the infinite volume 
limit. 

In the following two sections, we will discuss the repre- 
sentation for the repulsive and attractive models. Each 
will require a separate set of contours of integrations in 
order for the representation to be valid. We will notice 
that in the repulsive case, it is sufficient to integrate over 
the real line. The attractive case will require the use of 
contours separated out in the imaginary direction (to be 
qualified below) consistent with the fact that the spec- 
trum consists of "strings" with momenta taking values 
as complex conjugate pairs jig . 



A. Repulsive case 

We begin by discussing the repulsive case, c > 0. For 
this case, a similar approach has been independently de- 
veloped in Ref. [SJ and has been used by Lamacraft [33| 
to calculate noise correlations in the repulsive model. We 
will start with a generic initial state given by 

l*o)= Usix)l[bHx,m. (IV.2) 



with $s symmetrized. Using the symmetry of the bosonic 
operators, we can rewrite this state in terms of N-boson 
coordinate basis states, 



|*o) = N\ / $,(f)|x) 



(IV.3) 



where, 



the integral signs to obtain, 



\x)^e{x)Y[b\x,)\o). 



(IV.4) 



with 9{x) = 9{x\ > X2 > ■ ■ ■ > xn). It suffices therefore 
to show that we can express any coordinate basis state 
as an integral over the Bethe Ansatz eigenstates 

\x)^eix)ll[^Ai\,x)\X) (IV.5) 



\x,t)^9{x) I J]^ e— (^^)*A(A,f)|A). (IV.9) 



Attractive case 



with appropriately chosen contours of integration {7^} 
and A{X,x), which plays a role similar to the overlap of 
the eigenstates and the initial state. 

We claim that in the repulsive case equation (|IV.5|) is 
realized with 



AiX,x)=l[e- 






(IV.6) 



and the contours jj running along the real axis from 
minus to plus infinity. In other words eqn. pV.5[) takes 
the form, 






'y 



X 



J|e'^^(^^-"^)6^(yj)|0). (IV.7) 



Equivalently, we claim that the A integration above pro- 
duces lljS{yj - x-j). 

We shall prove this in two stages. Consider first y^ > 
Xn- To carry out the integral using the residue theorem, 
we have to close the integration contour in Ajv in the 



upper half plane. The poles in Aat are at \j 



i < N. 



These are all below j^ and so the result is zero. This 
implies that any non-zero contribution comes from y^ < 
Xn- Let us now consider yN-i > xn-i- The only pole 
above the contour is A^_j^ — Xn + *c. However, we also 
have, yN-i > xn-i > xn > yn =» yn-i > yn- 
This causes the only contributing pole to get canceled. 
The integral is again zero unless j/tv-i < xn^i- We can 
proceed is this fashion for the remaining variables thus 
showing that the integral is non-zero only for yj < Xj. 

Now consider j/i < xi. We have to close the contour 
for Ai below. There are no poles in that region, and the 
residue is zero. Thus the integral is non-zero for only 
2/1 = xi. Consider 7/2 < X2- The only pole below, at 
A2 = Ai — ic is canceled as before since we have 2/2 < 
X2 < xi ~ yi- Again we get that the integral is only 
non-zero for 2/2 ~ X2- Carrying this on, we end up with 



(IV.8) 



x) ll[S{y,^x,)bHy,)\0} 



In order to time evolve this state, we can act on it 
with the unitary time evolution operator. Since the in- 
tegrals are well-defined we can move the operator inside 



We now consider an attractive interaction, c < 0. 
As mentioned earlier, this changes the spectrum of the 
Hamiltonian, allowing complex, so-called string solu- 
tions, which in this model, correspond to many-body 
bound states. In fact, the ground state at T = consists 
of one iV-particle bound state. We will see how the Yud- 
son integral representation takes this into account. Sim- 
ilar properties are seen to emerge in [35|, where the au- 
thors obtain a propagator for the attractive Lieb-Liniger 
model by analytically continuing the results obtained by 
Tracy and Widom [33] for the repulsive model. 

One will immediately notice from the eigenstates that 
the pole structure of the S-matrix is altered. This change 
prevents the proof of the previous section from working. 
In particular, the poles in the variable A^v are at A-,- -l-?|c| 
for j < N, and the residue of the contour closed in the 
upper half plane is not zero any more. We need choose 
a contour to avoid this pole. This can be achieved by 
separating the contours in the imaginary direction such 
that adjacent Im[Aj — Xj-i] > \c\. At first sight, this 
seems to pose a problem as the quadratic term in expo- 
nent diverges at large positive A and positive imaginary 
part. There are two ways around this. We can tilt the 
contours as shown in Fig.[l]so that they lie in the conver- 
gent region of the Gaussian integral. The pieces towards 
the end, that join the real axes though essential for the 
proof to work at < = 0, where we evaluate the integrals 
using the residue theorem, do not contribute at finite 
time as the integrand vanishes on them as they are taken 
to infinity. Another more natural means of doing this 
is to use the finite spatial support of the initial state. 
The overlaps of the eigenstates with the initial state ef- 
fectively restricts the support for the A integrals, making 
them convergent. 



The proof of equation pV.5[) now proceeds as in the 
repulsive case. We start by assuming that yN > xn 
requiring us to close the contour in Aat in the upper 
half plane. This encloses no poles due to the choice 
of contours and the integral is zero unless yN < xn- 
Now assume yN-i > xn-i- Closing the contour above 
encloses one pole at A 



Af-^l 



A 



N 



i\c\, however since 
2/jv-i > Xn-i > Xn > yNi this pole is canceled by the 
numerator and again we have yN-i < xn-i- We proceed 
in this fashion and then backwards to show that the in- 
tegral is non-zero only when all the poles cancel, giving 
^s Ylj HVj ~ ^j)j as required. 




FIG. 4. Contours for the A integration. Shown here are three 
contours, and the closing of the TVth contour as discussed in 
the proof. 



V. TWO PARTICLE DYNAMICS 

We begin with a detailed discussion of the quench dy- 
namics of two bosons. As we saw, it is convenient to ex- 
press any initial state in terms of an ordered coordinate 
basis, \x) = 9{xi > X2 > ■ ■ ■ > XN)Yljb'^ixj)\0). At fi- 
nite time, the wave function of bosons initially localized 
at xi and X2 and subsequently evolved by a repulsive 
Lieb-Liniger Hamiltonian is given by. 
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where a = 2ct — i{yi — xi) — i(j/2 ~ X2). The above ex- 
pression retains the Bethe form of wave functions defined 
in different configuration sectors. The only scales in the 
problem are the interaction strength c and the scale from 
the initial condition, the separation between the particles 
at t = 0. 

In order to get physically meaningful results we need 
to start from a physical initial state. We choose first the 
state |\l!'(0)iatt) where bosons are trapped in a periodic 
trap forming initially a lattice- like state (see fig.[T^), 
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If we assume that the wave functions of neighboring 

bosons do not overlap significantly, i.e., e~°^ <^ 1, then 
the ordering of the initial particles needed for the Yudson 
representation is induced by the non-overlapping support 
and it becomes possible to carry out the integral analyt- 
ically. 

We now calculate the evolution of some observable 
in the state |\l/(0)iatt)- Consider first the evolution of 
the density p{x) = b^{x)b{x) at a; = 0. Fig. [5] shows 
(^iatt,i|p(0)|^iatt,^) for repulsive, attractive and non- 
interacting bosons. No difference is discernible between 
the three cases. The reason is obvious: the local inter- 
action is operative only when the wave functions of the 
particles overlap. As we have taken a <S^ a this will occur 
only after a long time when the wave-function is spread 
out and overlap is negligible. 
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FIG. 5. Color on line. {p{x = 0,i)) vs. t, after the quench 
from I'l/iatt). cr/a ~ 0.1. The curves appear indistinguishable 
since the particles start out with non significant overlap. The 
interaction effects would show up only when they propagated 
long enough to have spread sufficiently to reach a significant 
overlap, at which time the density is too low. 



Consider now an initial state where we set the sepa- 
ration a to zero, starting with maximal initial overlap 
between the bosons |vl/(0)cond)- We refer to this state as 
a condensate (in position-space) . Fig. [5] shows the den- 
sity evolution for attractive, repulsive and no interaction. 
The decay of the density is slower for attractive model 
than the for the non-interacting which in turn is slower 
than for the repulsive model - indeed, unlike before, the 
interaction is operative from the beginning. Still, the 
density does not show much difference between repul- 
sive and attractive interactions in this case. However, a 
drastic difference will appear when we study the noise 
correlations {'iQ,t\p{xi)p{x2)\'^o,t), as will be shown be- 
low. 

A comment about the attractive case is in order here. 
Recall that the contours of integration are separated in 
the imaginary direction. In order to carry out the inte- 
gration over A, we shift the contour for A2 to the real 
axis, and add the residue of the pole at A2 = Ai -I- i\c\. 
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FIG. 6. Color on line. (p(x = 0,f)) vs. t, after the quench 
from l^cond)- o" ~ 0.5, a = 0. As the bosons overlap interac- 
tion effects show up immediately. 



The two particle finite time state can be written as 
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7c refers to contours that are separated in imaginary 
direction, 7^ refers to all A integrated along real axis. 
/(A2 = Ai + i|c|) is the the residue obtained by shifting 
the A2 contour to the real axis from the pole at A2 at 
Ai +i|c|. This second term corresponds to a two-particle 
bound state. It is given by 

/(A2 = Ai +i|c|) = -2c I I e*^i('^i-^i)+'(^^-^=^) 

piAl (1/1-3:1 -1/2+^2) + §(y2-yi)4-f(a;i-X2) 
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where the last step involves symmetrizing the integrand 

in j/i and y2- Since c is negative, this term is a bound 

— Ml _ I 
state of the form e 2 W^ v^\ , Such bound states appear 

for any number of particles involved. For instance, for 
three particles, the Yudson representation with complex 
A's automatically produces multiple bound-states coming 
from the poles, i.e. /(A2 = Ai+i|c|),/(A3 = Ai-|-z|c|), etc. 
They give rise to two and three particle bound states of 
the form g-T^d^i-^'^l+l^'i-s'^l+l'ys-yal), The contribution 
to the above integral coming from the real axis corre- 
sponds to the non-bound state. 

Finally, putting together (|V.3|) and (|V.4[) . we get a 
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where a = —2\c\t — i(yi — xi) — i(j/2 — X2)- Surprisingly, 
the wave function maintains its form and we only need to 
replace c -^ —c. This simple result is not valid for more 
than two particles. 

We now compute the noise correlation in the evolving 
state. We expect the interaction to have a a significant 
effect as the geometry of the set up measures the inter- 
ference of "direct" and "crossed" measurements as shown 
in fig. [7^. In contrast, the density measurements do not 
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FIG. 7. (a) Color on line. The Hanbury-Brown Twiss effect, 
where two detectors are used to measure the interference of 
the direct and the crossed waves. The S'-matrix enters explic- 
itly, (b) The density measurement is not directly sensitive to 
the S-matrix. The thick black line shows the wave-function 
amplitude, the dotted lines show time propagation. 

see the S'-matrix, as shown in fig. [7f). 

This is the famous Hanbury-Brown Twiss experi- 
ment [37[ where for free bosons or fcrmions, the cross- 
ing produces a phase of ±1 and causes destructive or 
constructive interference. In our case the set up is gener- 
alized to multiple time dependent sources with the phase 
given by the two particle S'-matrix capturing the interac- 
tions between the particles. In Fig.|5]we present the two 
point correlation matrix {p{xi)p{x2)) for the repulsive 



gas, attractive gas, and the non-interacting gas, shown 

at different times, starting with the lattice initial state. 

Figure [9] shows the same for the condensate initial state. 

In both these, we note that the repulsive gas develops 
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FIG. 9. Color on line. Time evolution of density-density 
correlation matrix {{p{x)p{y))) for the l^&cond) initial state. 
Blue is zero and red is positive. The repulsive model shows 
anti-bunching, i.e., fermionization at long times, while the 
attractive model shows bunching. 



FIG. 8. Color on line. Time evolution of density-density 
correlation matrix {{p{x)p{y))) for the |^iatt) initial state. 
Blue is zero and red is positive. The repulsive model shows 
anti-bunching, i.e., fermionization at long times, while the 
attractive model shows bunching. 

fcrmionic correlations (i.e., strong anti-bunching), and 
the attractive gas retains bosonic correlations at long 
time, showing strong bunching. 

It is interesting to compare this result with the time 
evolution after a quench on the lattice by the Bose- 
Hubbard model, the lattice counterpart of the Licb- 
Linigcr model, as we shall see in Appendix |Al The re- 
sults for continuum model differs strongly from those of 
the lattice model. 

We expect the results to be qualitatively similar for 
higher particle number. In order to go beyond two parti- 
cles however, the integrations cannot be carried out ex- 
actly. However, we can extract the asymptotic behavior 
of the wavefunction analytically, as we show below. 



inant. The other regime where N is sent to infinity first 
will be discussed in a separate report. We first deal with 
the repulsive model, for which no bound states exist and 
the momentum integrations can be carried out over the 
real line, and then proceed to the attractive model. In 
a separate sub-section, we examine the effect of starting 
with a condensate-like initial state. 



A. Repulsive interactions - Asymptotics 



VI. MULTIPARTICLE DYNAMICS AT LONG 
TIMES 

In this section, we derive an expression for the multi- 
particle wavefunction evolution at long times. The num- 
ber of particles N is kept fixed in the limiting process, 
hence, as discussed in the introduction we are in the low 
density limit where interactions arc expected to be dom- 



From (|IV.5[) we can sec by scaling A -^ A-^t, we get 



^f,(A. - A,) ^ sgn(y, _ y^) + O ( _ ) , (VI.l) 
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yielding to leading order, 
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c^(y) being fermionic creation operators replacing the 
"fermionized" hardcore bosonic operators, Yij'^Hyj) = 

Ylt<j sgn(yj - yj)b^yj). We denote h/^ = J^ dc^{x)dc{x) 
the free fermionic Hamihonian and Ay is an anti- 
symmetrizer acting on the y variables. Thus, the re- 
pulsive Bose gas, for any value of c > 0, is governed 
in the long time by the c = oo hard core boson limit 
(or its fermionic equivalent) [3l|, [Sg, and the system 
equilibrates with an asymptotic momentum distribution, 
"A: = (^o|c|,Cfe|^o), determined by the antisymmetric 
wavefunction ^oiv) = Ay9{y)'i>o{y) and the total energy. 

We will now derive the corrections to the infinite time 
limit. At large time, we use the stationary phase approxi- 
mation to carry out the A integrations. The phase oscilla- 
tions come primarily from the exponent e~'^j t+t^j (vj -a^j ) 
At large t (i.e., t ^ -ij), the oscillations are rapid, and 
the stationary point is obtained by solving 



d 
dX 



[—iX^t + iXjivj — Xj)] = 0. 



(VI.2) 



Note that typically one would ignore the second term 
above since it doesn't oscillate faster with increasing t, 
but here we cannot since the integral over y produces 
a non-zero contribution for j/ ~ t at large time. Doing 
the Gaussian integral around this point (and fixing the 
5- matrix prefactor to its stationary value) , we obtain for 
the repulsive case. 
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In the above expression, the wavefunction has support 
mainly from regions where yj/t is of order one. In an 
experimental setup, one typically starts with a local finite 
density gas, i.e., a finite number of particles localized over 
a finite length. With this condition, at long time, we can 



neglect Xj/t in comparison with yj/t, giving 
\x,t) -^ fY[z,,i^,-i,)Y[^/'^'^-^^^^^bHy,)\0) 
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(VI.4) 
where i, ~ jz- 

We turn now to calculated the asymptotic evolution 
of some observables. To compute the expectation value 
of the density we start from the coordinate basis states, 
{x' ,t\p{z)\x,t) which we then integrate with the chosen 
initial state, 

{x',t\p{z)\x,t)^Y. f \Y.^(y^-'n 
{p} •'y \ J ) 



i<j 



(VI.5) 



Note that the above product of iS-matrices is actually 
independent of the ordering of the y. First, only those 
terms appear in the product for which the permutation 
P has an inversion. For example, say for three particles, 
if P = 312, then the inversions are 13 and 23. It is only 
these terms which give a non-trivial iS-matrix contribu- 
tion. For the non-inverted terms, here 12, we get 

gi ^ 6 - ic8gn{yi - y2) Ci - 6 + ^csgn(yi - t/2) 

^1 - 6 - «c ^1 - 6 + ic 

(VI.6) 
which is always unity irrespective of the ordering of yi , 1/2 ■ 
For a term with an inversion, say 23, we get, 

6 - 6 - icsgn(y2 - ys) 6 - 6 + icsgn(?/3 - y2) 



6 - 6 - ic 
which is always equal to 

6 - 6 - ic 



6 - 6 + ic 
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irrespective of the sign of 2/2 ~ 2/3- This allows us to carry 
out the integration over the yj. 



1. Lattice initial state 

In order to calculate physical observables, we have 
to choose initial states. We choose two different initial 
states for the problem, one with N particles distributed 
with uniform density in a series of harmonic traps given 
by, 

f JL 1 (xj+(,-l)o)^ 

l^iatt) = / n^^TT^ ^ bHx^m, (VI.9) 

such that the overlap between the wave functions of two 
neighboring particles is negligible. In this particular case. 
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the ordering of the particles is induced by the hmited 
non-overlapping support of the wave function. 

In the lattice-like state, the initial wave function starts 
out with the neighboring particles having negligible over- 
lap. At small time (as seen from (|V.1[) ). the particle re- 
pel each other, but they never cross due to the repulsive 
interaction. So at large time, the interaction docs not 
play a role since the wave functions are sufficiently non- 
overlapping. It is only the P = 1 contribution then that 
survives, and we get for the density 
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We need to integrate the position basis vectors \x) over 
some initial condition. We do this here for the lattice 
state (|VI.9P This gives 



(*latt,i|p(^)|^latt,i) = /Olatt(6) 
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Mathematically, any S'-matrix factor that appears will 
necessarily have zero contribution from the pole - this is 
easy to see from the pole structure, and the ordering of 
the coordinates. In order to get a non-zero result, we 
need to fix at least two integration variables (i.e., the 
yj). Thus the first non-trivial contribution comes from 
the two-point correlation function. 

We now proceed to calculate the evolution of the noise, 
i.e., the two body correlation function piiz , z' \t\:^tt = 
(^iatt,i|p(^)p(^')l^iatt,i)- The contributions can be 
grouped in terms of number of crossings, which corre- 
sponds to a grouping in terms of the coefficient e~'^° [34| . 
The leading order term can be explicitly evaluated and 
we show below which terms contribute. In general we 
have 

P2 iatt(z, z';t) = ^ ^ S{yj - z)5{yu - z') 

{P}-'y \3,K J 

X n SiC,-Q'[[^^'^^^''^-^-^^^\ (VI.12) 
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The above shorthand in the S'-matrix product means that 
only the (ij) that belong to the inversions in P are in- 
cluded. We will now determine which terms contribute in 
this sum. First note that for integration over a particular 
^j, the residue depends on the sign of Xj — xp.. Let us 
consider a specific example. Consider the three particle 
case with the term P — 321. All three S'-matrix factors 
appear in this term. ^3 has a pole at ^1 — ic and ^2 ~ ic. 
Thus integrating over j/3 will give a non-zero residue only 
if X3 > x\ which is however not satisfied by the initial 
conditions we choose. So, everything is zero, unless we 
do not integrate over 2/3, implying it has to be one of 
the measured variables. Similarly for ^i, the poles are at 
^2 + *c and ^3 -|- ic. In order to get a non-zero residue we 
need x\ < Xg which again is not satisfied by the initial 



conditions. We get a non-zero result if we pin ^1. As for 
^2 it has poles both above and below the real line and so 
this always gives a non-zero contribution irrespective of 
the sign of X2 — x^- 

One can see that this argument can be extended to 
the case with more particles. Depending on what coordi- 
nates we are measuring at, we'll get a specific contribu- 
tion from the sum over permutations. The next simplifi- 
cation comes from not allowing any crossings among the 
unmeasured coordinates. It can be shown that allowing 
for these gives us a higher order contribution in e^'^". In 
other words, the leading order contribution comes from 
terms such as P = 21, 32, 321, 4231, 5342, 52341, . . .. The 
only exchanges are on the ends. A general term will 
therefore look like (for I < k), 

I 5{yi~z)5{yu^z') \{ S(e. - 0)^(0 - 6) 
•'y 3^{i.k} 
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We have to sum over Z, fc, which will automatically sum 
over the number of intermediate j's appearing. We'll 
integrate the above general term, since the yj integrals 
factor anyway. This gives 
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We can sum the different contributions now. Note that 
the number of terms appearing the product over j is given 
by fc — ^ — 1 . So for a given I , we have to sum over all the 
k. Using a shorthand notation, the sum can be written as 
(note that it is understood that yi and yk are integrated 
over using the delta functions. We retain the indices to 
keep track of the terms. We actually have £1 = 2zt and 
6 = 2z'i. 
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In order to proceed with the summation, we have to inte- 
grate over the x. We use the initial conditions described 
by (jVI.Qp . i.e., the lattice-like state. We'll do it term by 
term. 
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Note that last expression has no j dependence. So, the 
sum over the product over j is just a geometric series. 
Recah that ^i and ^fc arc fixed at z and z' respectively. 
This series can be summed. Note that the number of 
terms in the product is equal to the k — I, and summing 
over them is effectively summing over k. The previous 
term therefore needs to be taken into account. Writing 
6 = Cz and ^k — ^z' , we can write the sum as 
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Finally, we have to account for the k > / case which is 
equivalent to setting ^j = ^^i and ^i = £,z- Doing this is 
further equivalent to adding the complex conjugate. We 
also have to take into account the term with no permu- 
tations. So, finally we have. 
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To compare with the Hanbury-Brown Twiss result, we 
calculate the normalized spatial noise correlations, given 
by C2{z,z') = J0£l _ 1 ^ C2{z,z'). In the non- 

interacting case, i.e., c = 0, S{£) = 1 and gzz' = and 
we recover the HBT result for A^ = 2, 



C2°(C.,e.') = ^cos(a(e.-e;)) (VI.22) 

One can also check that the limit of c — )• 00 gives the 
expected answer for free fermions, namely. 



CTi^z.iz-) 



1 



cos(a(^;, - Q) 



(VI.23) 



At finite c we can see a sharp fermionic character appear 
that broadens with increasing c as shown in Fig. 1101 The 
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FIG. 10. Color online. Normalized noise correlation function 
C2{£,; — C)- Fermionic correlations develop on a time scale r 
c 

Ref.ll 



, so that for any c we get a sharp fermionic peak near 
~ 0, i.e., at large time. The key shows values of ca (from 



large time behavior is captured in a small window around 
^ = 0. One can see that at any finite c, the region near 
zero develops a strong fermionic character, thus indicat- 
ing that irrespective of the value of the coupling that we 
start with, the model flows towards an infinitely repul- 
sive model at large time, that can be described in terms 
of free fermions. We also obtained this result "at" t = 00 
at the beginning of this section. 

For higher particle number, we see "interference 
fringes" corresponding to the number of particles, that 
get narrower and more numerous with an increase, mem- 
ory of the initial lattice state. However, the asymptotic 
fermionic character does not disappear. Figures [TT] and 
1121 show the noise correlation function for five and ten 
particles respectively. The large peaks are interspersed 
by smaller peaks and so on. This reflects the character 
of the initial state. 




FIG. 11. Normalized noise correlation function for five parti- 
cles released for a Mott-like state for c > (from Ref. llTl ) 



2. Quenching from a bound state 

In this brief section our initial state is the ground state 
of the attractive Lieb-Liniger Hamiltonian (with interac- 
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FIG. 12. Normalized noise correlation function for ten parti- 
cles released for a Mott-like state for c > 0. 



tion strength — cq < 0. For two bosons, this take the 
form \M, 
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(VI.24) 
and we quench it with a repulsive Hamiltonian.Thc long 
time noise correlations are displayed in Fig. [T31 We see 
that while the initial state correlations are preserved over 
most of the evolution, in the asymptotic long time limit 
the characteristic fermionic dip. We expect similar effects 




FIG. 13. Normalized noise correlation function for two par- 
ticle quenched from a bound state into the repulsive regime. 
The legend indicates the values of c that the state is quenched 
into. We start with cqct^ =3, a = 1, cq being the interaction 
strength of the initial state Hamiltonian. Again, we see the 
fermionic dip, but the rest of the structure is determined by 
the initial state. 



them, leading to sum over several terms. Fig. 1141 shows an 
example of how this works. 
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FIG. 14. Contribution from stationary phase and pole at large 
time in the attractive model. The blue curve represents the 
shifted contour. 

In Ref. llTl a formula was provided for the asymptotic 
state. Here we give a more careful treatment by taking 
into account that the fixed point of the approximation 
moves for terms that come from a pole of the ^-matrix. 
It is therefore necessary to first shift the contours of in- 
tegration, and then carry out the integral at long time. 
We carry this out below. 

Shifting a contour over a pole leads to an additional 
term from the residue: 



f dX2 f dA2 .^,. 

•'12 -"" -^72 



Al+^|c|) (VI.25) 



where TZ{x) indicates that we evaluate the residue given 
by the pole x. jj indicates the original contour of inte- 
gration and 7^ indicates that integration is carried out 
over the real axis. Proceeding with the other variables 
we end up with 



•^7li72,--- ,7jv Jjf 



+in{X2 -^ Ai +i\c\[ 
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+i7^(A3 ^ Ai + i\c\) + in{\3 ^ Aa + i\c\) 



73 



7S 



+in{XN -^ Xi+ i\c\) + in{XN ^ A2 + i\c\)+ 

-|-•••+i7^(AAr ^ AAr_i+i|c|)] (VI.26) 



The integrals can now be evaluated using the stationary 
phase approximation. The correction produced by the 
above procedure does not affect the qualitative features 
observed in Ref. [13. 



for any number of bosons. 



Lattice initial state 



B. Attractive interactions 

For the attractive case, since the contours of integra- 
tion are spread out in the imaginary direction, we have 
the contributions from the poles in addition to the sta- 
tionary phase contributions at large time. The stationary 
phase contribution is picked up on the real line, but as we 
move the contour, it stays pinned above the poles and we 
need to include the residue obtained from going around 



We now calculate the evolution of the density and the 
two body correlation function in order to compare with 
the repulsive case. We will first study the two particle 
case. Although we have a finite time expression for this 
case from which we can directly take a long time limit, we 
will study the asymptotics using the above scheme for an 
A^-particle state, since we have an analytical expression 
to go with. We get two terms, the first being the sta- 
tionary phase contribution, and is just like the repulsive 
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case with c — > — c. The second is the contribution from 
the pole. It contains the bound state contribution which 
brings about another interesting feature of the attractive 
case. While the asymptotic dynamics of the repulsive 
model is solely dictated by the new variables ^j = ||, 
and all the time dependence of the wave function enters 
through this "velocity" variable, this is not the case in 
the attractive model. While it is true that the system 
is naturally described in terms of ^ variables, there still 
exists non-trivial time dependence. 

First, we integrate out the x dependence assuming an 
initial lattice-like state. This gives. 




I* 






y/Anit 
(VI.27) 



FIG. 15. Color on line. Variation of C2 for the attractive case 
with time. Note the growth of the central peak. 



Defining ^i^t) from \^uttit)) = /^ •/'(C, 11, ^^(%)|0), 
we have for the density evolution under attractive inter- 
actions, c < 0, 



Pla 



(z;i) = E [S(yj - zW{^P,m{U) (VI.28) 



{P}j 



We can show numerically (the expressions are a bit un- 
wieldy to write here), that asymptotically, the density 
shows the same Gaussian profile that we expect from a 
uniformly diffusing gas, namely, e~* '^ . 

With this, we can proceed to compute the noise cor- 
relation function. The two particle case is easy, as there 
are no more integrations to carry out. We get. 
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FIG. 16. Color on line. C2{£,,~£,) for three particles in the 
attractive case plotted for three different times. At larger 
times, the correlations away from zero fall off. ta^ = 20, 40, 60 
z)S{yk - z')(t>* (^p, t)4>{^, t) for blue, magenta and yellow respectively, (from Ref. [l7h 



(VI.29) 

where ^g is the symmetrized wavefunction. FiglTSJ shows 
the normalized noise correlations for different values of t 
and Figure [!§] shows the variation with c. 

For more particles, we see interference fringes similar 
to the repulsive case. We note that the central peak 
is now increasing and sharpening with time, indicating 
increasing contribution from bound states (see Fig. [12] 
for an example). 



C. Starting with a condensate - attractive and 
repulsive interactions 

In this section, we study the evolution of the Bose gas 
after a quench from an initial state where all the bosons 
are in a single level of a harmonic trap. For t < 0, the 
state is described by 



Recall that in order to use the Yudson representation, 
the initial state needs to be ordered. We can rewrite the 
above state as 



I*. 



(xi > 



> xn) 



s^U 



{ita^) 



rbHx.m (VI.31) 



where 5 is a symmetrizer. The time evolution can be 
carried out via the Yudson representation, and again, we 
concentrate on the asymptotics. For the repulsive model, 
the stationary phase contribution is all that appears, and 
we get 



n^ 



y t<j 



Hi yj 3Ci -r Xj 
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(VI.30) 
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At large time t, we therefore have 



I*. 



i(i)> 



(x'l > • • • > XN)(f>2{x)I{y, X, t) 

xllb\y,)\0) (VI.33) 



4>2{x) is symmetric in x. I{y,x,t) is symmetric in the y 
but not in the x. Therefore we have to carry out the x 
integration over the wedge xi > ■ ■ ■ > xm- This is not 
straightforward to carry out. If I{y,x,t) was also sym- 
metric in X, then we can add the other wedges to rebuild 
the full space in x. However, due to the S'-matrix factors, 
symmetrizing in y does not automatically symmetrize in 
X. The exponential factors on the other hand are au- 
tomatically symmetric in both variables if one of them 
is symmetrized because their functional dependence is of 



the form f{yj 



It is however possible to make the S'- 



matrix factors approximately symmetric in x, and we will 
define what we mean by approximately shortly. What is 
important is to obtain a yj — Xj dependence. As of now, 
the S'-matrix that appears in the above expression is 






Vi - Vj 



2t 



2t 



icsgniy, - y^) 



Vi-Vj 



-Vi-Xi+Xi 



(VIM) 

First, we can change sgii{yi — yj) to sgn ( ^'~^^ j since t > 

0. Next, note that asymptotically in time, the station- 
ary phase contribution comes from ^ ~ C'(l). However, 
since x has finite extent, at large enough time, ^ 

We are therefore justified in writing sgn I ^ 
The only problem could arise when yi 



0. 



Vj- 



2t 



2t 

yj . How- 
ever, if this occurs, then the S-matrix is approximately 
sgn(yi — yj) which is antisymmetric in ij. With this pref- 
actor the particles are effectively fermions, and therefore 
at yi ~ yj, the wave- function has an approximate node. 
At large time therefore, we do not have to be concerned 
with the possibility of particles overlapping, and includ- 
ing the Xi inside the sgn function is valid. With this 
change the S-matrix also becomes a function of yj — Xj 
and symmetrizing over y one automatically symmetrizes 
over X. 

In short, we have established that the wave function 
asymptotically in time can be made symmetric in x. This 
allows us to rebuild the full space. We get 



I*. 



x,V p 



p j 

M^)r{y,x,t)Y[b^iy,)\0) 

3 

(VI.35) 



where the s superscript indicates that we have established 
that /(y, X, t) is also symmetric in x. With this in mind, 
we can do away with the ordering when we're integrating 
over the x if we symmetrize the initial state wave function 



and the final wave function. Note that when we calcu- 
late the expectation value of a physical observable, the 
symmetry of the wavefunction is automatically enforced, 
and thus taken care of automatically. 

Recall that when we calculated the noise correlations 
of the repulsive gas, in order to get an analytic expression 
for N particles, we considered the leading order term, i.e., 
the HBT term. We did this by showing that higher or- 
der crossings produced terms higher order in e"^'^" which 
we claimed was a small number. Now, however, a = 0, 
and although the calculation is essentially the same with 
our approximate symmetrization, this simplification does 
not occur. The two and three particle results remain an- 
alytically calculable, but for higher numbers, we have to 
resort to numerical integration. Fig. [T7] shows the noise 
correlation for two and three repulsive bosons starting 
from a condensate. For non-interacting particles, we ex- 
pect a straight line C2 = h- When repulsive interactions 
are turned on, we see the characteristic fermionic dip de- 
velop. The plots for the attractive Bose gas are shown 
in Figs. [TSl and [T9I As expected from the non- interacting 
case the oscillations arising from the interference of par- 
ticles separated spatially does not appear. The attractive 
however docs show the oscillations near the central peak 
that are also visible in the case when we start from a 
lattice-like state. It is interesting to note that for three 
particles we do not see any additional structure develop 
in the attractive case. 




FIG. 17. C2(C, —0 fo'" l"*^o (blue) and three repulsive bosons 
starting from a condensate. Unlike the attractive case, there 
is no explicit time dependence asymptotically, ca = 3 (from 
Ref . 113) 



VII. CONCLUSIONS AND THE DYNAMIC RG 
HYPOTHESIS 

We have shown that the Yudson contour integral rep- 
resentation for arbitrary states can indeed be used to 
understand aspects of the quench dynamics of the Lieb- 
Linigcr model, and obtain the asymptotic wave functions 
exactly. The representation overcomes some of the major 
difficulties involved in using the Bethe-Ansatz to study 
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FIG. 18. Noise correlation for two attractive bosons start- 
ing from a condensate - as time increases, the central peak 
dominates 
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FIG. 19. C2(C, — for three attractive bosons starting from 
a condensate. Note that the side peak structure found in 
fig. [16] is missing due to the initial condition. We show the 
evolution at three times. As time increases, the oscillations 
near the central peak die out. Times from top to bottom 



tc^ = 20, 40, 60. (from Ref. [17 



the dynamics of some integrable system.s by automati- 
cally accounting for complicated states in the spectrum. 

We sec some interesting dynamical effects at long 
times. The infinite time limit of the repulsive model 
corresponds to particles evolving with a free fermionic 
Hamiltonian. It retains, however, memory of the initial 
state and therefore is not a thermal state. The correlation 
functions approach that of hard core bosons at long time 
indicating a dynamical increase in interaction strength. 
The attractive model also shows a dynamic strengthening 
of the interaction and the long time limit is dominated 
by a multiparticle bound state. This of course does not 
mean that it condenses. In fact the state diffuses over 
time, but remains strongly correlated. 

We may interpret our results in terms of a "dynamic 
RG" in time. The asymptotic evolutions of the model 
both for c > and for c < are given by the Hamil- 
tonians H^ with c — >■ ±oo respectively. Accepting the 
RG logic behind the conjecture one would expect that 
there would be basins of attraction around the Lieb- 



Liniger Hamiltonian with models whose long time evolu- 
tion would bring them close to the " dynamic fixed points" 
Hj.. One such Hamiltonian would have short range po- 
tentials replacing the (5-function interaction that renders 
the Lieb-Liniger model integrable. Perhaps, lattice mod- 
els could be also found in this basin whose time asymp- 
totics would be close, in the repulsive case, to that given 
by a free fermionic model on the lattice. Clearly, as dis- 
cussed earlier, the Bose-Hubbard model is not such a 
model since it has a lattice symmetry that is not present 
in the Lieb-Liniger model. This could be however over- 
come by adding such terms as the next nearest hopping 
or interactions that break this symmetry, or as shown in 
Appendix |X1 with an appropriate choice of initial state. 

We have to emphasize, however, that as these models 
arc not integrable, we do not expect that they would ac- 
tually flow to H^. Instead, starting close enough in the 
"basin", they would flow close to H^ and spend much 
time in its neighborhood, eventually evolving into an- 
other, thermal state. We thus conjecture that away from 
integrability, a system would approach the corresponding 
non-thermal equilibrium, where the dynamics will slow 
down leading to a "prethermal" state [33 • Fig. [^ shows 
a schematic of this. Such prethermalization behavior has 
indeed been observed in lattice models [40|- The system 
is expected to eventually find a thermal state. It is there- 
fore of interest to characterize different ways of breaking 
integrability to see when a system is "too far" from inte- 
grability to see this effect and in what regimes a system 
can be considered as close to integrabilit y. F or a review 
and background on this subject, see Ref. |4l|. 

Further, the flow diagram in Fig. [201 might have an- 
other axis that represents initial states. Studying the 
Bose-Hubbard model shows an interesting initial state 
dependence. Whereas the sign of the interaction does 
not affect the quench dynamics, the asymptotic state de- 
pends strongly on the initial state, with a lattice-like 
state leading to fermionization, and a condensate-like 
state retaining bosonic correlations. The strong depen- 
dence on the initial state in the quench dynamics is ev- 
ident from eq. (jl.2[) and is subject of much debate, in 
particular as relating to the Eigenstate Thermalization 
Hypothesis [l,!!!!!!]. 
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FIG. 20. Schematic showing pre-thermalization of states in a 
non-integrable model 

This work also opens up several new questions. It 
provides a prediction for experiments that can be car- 
ried out in the context of continuum cold atom systems 
(though the experiments we are aware of are carried 
out on the lattice and therefore described by the Bose- 
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Hubbard model) Theoretically, while the representation 
is provable mathematically, further investigation is re- 
quired to understand, physically, how it achieves the te- 
dious sum over eigenstates, while automatically account- 
ing for the details of the spectrum. This would allow us 
to extend the approach to other models with a more com- 
plicated S'-matrix structure. It would also be useful to 
tie this approach to other means of calculating overlaps 
in the Algebraic Bethe Ansatz, i.e., the form- factor ap- 
proach. The representation can essentially be thought of 
as a different way of writing the identity operator. From 
that standpoint, it could serve as a new way of evalu- 
ating correlation functions using the Bethe Ansatz. We 
are currently studying generalizations of this approach to 
other models that can be realized in optical lattices. 
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Appendix A: Quenching the Bose-Hubbard model 

We compare the results obtained in Section |V] with 
those from the lattice version of the Lieb-Liniger model 
- the Bose-Hubbard model. 



^BH — 



n 



tb]h+i + h.c. + Un,{n, - 1) (A.l) 



It describe bosons b hopping on a Id lattice with on- 
site interaction U and is non-integrable since it allows 
multiparticle interactions on the same site. It has been 
extensively studied in many contexts and much is known 
about its equilibrium properties (see e.g., Ref. |4J). For 

< U /t <^ 1, the model is a superfluid, and for U /t ^ 

1 it is a Mott insulator. For negative U , the model is 
attractive and the ground state is a Bose condensate. 
A non- equilibrium phase diagram of the Bose-Hubbard 
model is given in Ref. |45| . 

We study here the two boson quench dynamics and 
contrast it with the corresponding dynamics of the Lieb- 
Liniger model. Contrary to what one may expect, the 
introduction of the lattice modifies the dynamics in an 
essential way even at long times and distances. The cal- 
culations of density correlations as a function of time af- 
ter a sudden quench have been carried out using the Al- 
gorithms and Libraries for Physics Simulations (ALPS) 
code [H-Iiij and the Open source TEBD package [ii] af- 
ter making the necessary modifications to accommodate 
the initial states wc are interested in. Our results confirm 
some results obtained in Ref. [l^. 

In Fig. [21] we show the time evolution of the corre- 
lation matrix defined as {riirij) after a sudden quench 
from an initial state 6q&q|0), and in Fig. [^ the evolution 



from initial state ^JfoJlO). We quench into the interacting 
regime, where \U\/t ~ 10. There are a couple of interest- 
ing features: (1) Unlike the situation in the Lieb-Liniger 
model where the bunching or anti-bunching effect is inde- 
pendent of the initial state, here, quenching a lattice-like 
state leads to anti-bunching Fig. [211 while quenching a 
condensate-like state leads to bunching Fig. [5TJ It is also 
interesting to compare the the anti-bunching evolution of 
the bosons with the evolution of free fermions in Fig.[23l 




FIG. 21. Color online. Time evolution of the correlation 
matrix after a sudden quench from a state containing two 
bosons on the same site. The values increase from blue (0) to 
red. The correlations remain strong in the center indicating 
strong bunching. 

(2) The sign of the interaction plays no role in the evolu- 
tion in the Bose-Hubbard model, as seen from either fig- 
ures. This is unlike the situation in the continuum model 
where for repulsive interactions anti-bunching (fermion- 
ization) occurs independently of the initial state, while 
bunching will take place for attractive interactions. 

This non-dependence on the sign of the interaction is 
due to a particle-hole symmetry that is present on the 
lattice, but not in the continuum. The Id lattice, be- 
ing bipartite, allows the the transformation bj — )■ e^'^^bj, 
withj the site index, under which the hopping terms pick 
up a minus sign, while the on-site interaction terms are 
unaffected, thus U/t — > —U/t. In terms oiU the corre- 
sponding unitary operator we have. 



UHBHit,-U)U^ ^-HBnit^U). 



(A.2) 



Denoting eigenstates and eigenvalues of 7Jbh(^, C^) = 
Hbh by \m) and em and the corresponding eigenstates 
and eigenvalues of HBnit, —U) = -ffen by \m) and e™ we 
can relate the states by: \m) ~ U\m), and the eigenval- 
ues by: em = —^m- The time evolution of an operator O 
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FIG. 22. Color on line. Time evolution of the correlation 
matrix after a sudden quench from a state containing bosons 
on two neighboring sites. The values increase from blue (0) 
to red. The off diagonal correlations indicate anti-bunching, 
as can be seen from free fermion evolution in Fig. 1231 




FIG. 23. Time evolution of free fermions on a lattice. Notice 
how off diagonal correlations develop. The values increase 
from blue (0) to red. 



under the action of HBii{t, U) from an initial state |*I'o) 
{Oit))H = J2 (*o|m')(m|*o)(m'|0|m)e-'(^'"-^'"')* 

(A.3) 
Under Hbh, 

{0{t))^ = Y^ (*o|m')(m|*o>(m'|0|m)e-^('^--^"')* 

(A.4) 



Both initial state we considered, |^iatt) = ^o^i|0) and 
I'I'cond) = (&o)^|0) are simply transformed, U\"^o) = 
±|^o) and as they occur twice in the overlaps the trans- 
formation leaves no effect. Similarly the operators we 
have considered (density-density correlations) are bilin- 
ear in the site operators and are therefore not affected, 
U'^OU = 0. This gives 



(A.5) 
Next, we note that the Bosc-Hubbard Hamiltonian is in- 
variant under time reversal, i.e., [-ffBH,7^ = where T 
is the anti-unitary time reversal operator: 



TtT-^ = -t, TiT-' = -i. 



(A.6) 



Applying the time reversal operator to the expectation 
value above, we get 






m) e 



«(eT7i— e^/)t 



Y (*o|m)(7n'|*o)(m|0|m')e 



j(em-e,„')t 



m,rn' 



(A.7) 



= ;^(*o|m')(m|*o)(m'|0|TO)e 

771, m' 
= {0{t))H 



-2( €™ — €„ 



thus indeed, the time evolution looks the same for both 
signs of the interaction. Note that with initial states or 
operators that are not invariant (up to a sign) under the 
transformation U, we should see a difference in the time 
evolution of the attractive and repulsive models. 

A similar symmetry exists in the XXZ model or the 
Hubbard model in Id (or higher dimensional bipartite 
lattices). For the magnet, the sign of the anisotropy A 
leads to either ferromagnetic or antiferromagnetic ground 
states for negative or positive anisotropy. However, it 
does not influence the quench dynamics [19|, as can be 
seen from arguments like the above. Similarly in the 
Hubbard model, the quench dynamics is unaffected by 
the change of sign of U [50| • 
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